* Figure VII panels a)-b)
* 29 Oct 2021

clear all
set more off

grstyle init
grstyle set plain, nogrid horizontal
grstyle color background white

*Choose figure width:
local fig_w= 5
local res = `fig_w' * 600

import excel "$output_s/w3_test_of_model.xlsx", sheet("Sheet1") firstrow clear

* Bar graph of cases - model vs prediction
gen occup_opt = ceil(case_opt_w3/2) //Quick way to separate cases
gen occup = ceil(case_real_w3/2)

label define occupation 1 "Mixing" 2 "Livestock Rearing Only" 3 "Wage Work Only" 4 "No Work"
label values occup* occupation

drop if occup == 4

preserve
	keep hhid occup*
	rename occup_opt occup2
	rename occup occup1
	reshape long occup, i(hhid) j(tag)
	label define tags 1 "Actual" 2 "Predicted"
	label values tag tags
	catplot tag occup, percent(tag) var1opts(label(labsize(small))) var2opts(label(labsize(small))) recast(bar) asyvars bargap(20) bar(1, color(pink*0.4)) bar(2, color(green*0.4)) ytitle("Percentage") name(bars, replace) yla(0(10)30, gmax)
	graph save bars "$output/Bars_cases_excl_no_work_w3.gph", replace
	graph export "$output/Bars_cases_excl_no_work_w3.png", as(png) replace
restore


* Figure VII panel a)

	*Colored version:
	twoway (lpolyci l_opt_w3 k_w3 if k_w3<16000, lcolor(green) fcolor(green%30) alcolor(%30) bw(800) level(95)) (lpolyci l_w3 k_w3 if k_w3<16000, lcolor(pink) fcolor(pink%30) alcolor(%30) bw(800) level(95)), xtitle("Capital at Wave 3 (BDT)", margin(t+2)) ytitle("Hours in Livestock Rearing") yla(0(200)800) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_l, replace)
graph save Lpoly_l "$output/Lpoly_l_kw3.gph", replace
graph export "$output/Lpoly_l_kw3.png", as(png) replace
	
	*TIFF Greyscale version:
	twoway (lpolyci l_opt_w3 k_w3 if k_w3<16000, lcolor(gs10) lp(dash) fcolor(gs10%30) alcolor(%30) bw(800) level(95)) (lpolyci l_w3 k_w3 if k_w3<16000, lcolor(gs3) fcolor(gs3%30) alcolor(%30) bw(800) level(95)), xtitle("Capital at Wave 3 (BDT)", margin(t+2)) ytitle("Hours in Livestock Rearing") yla(0(200)800) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95%% CI")) name(Lpoly_l, replace) xsize(`fig_w')
	graph export "$output/TIFF/Lpoly_l_kw3.tif", as(tif) replace width(`res')

	
* Figure VII panel b)

	*Colored version:
	twoway (lpolyci h_opt_w3 k_w3 if k_w3<16000, lcolor(green) fcolor(green%30) alcolor(%30) bw(800) level(95)) (lpolyci h_w3 k_w3 if k_w3<16000, lcolor(pink) fcolor(pink%30) alcolor(%30) bw(800) level(95)), xtitle("Capital at Wave 3 (BDT)", margin(t+2)) ytitle("Hours in Wage Work") yla(0(200)1400) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_h, replace)
graph save Lpoly_h "$output/Lpoly_h_kw3.gph", replace
graph export "$output/Lpoly_h_kw3.png", as(png) replace

	*TIFF Greyscale version:
	twoway (lpolyci h_opt_w3 k_w3 if k_w3<16000, lcolor(gs10) lp(dash) fcolor(gs10%30) alcolor(%30) bw(800) level(95)) (lpolyci h_w3 k_w3 if k_w3<16000, lcolor(gs3) fcolor(gs3%30) alcolor(%30) bw(800) level(95)), xtitle("Capital at Wave 3 (BDT)", margin(t+2)) ytitle("Hours in Wage Work") yla(0(200)1400) legend(order(4 "Actual" 2 "Predicted" 3 "95% CI" 1 "95% CI")) name(Lpoly_h, replace) xsize(`fig_w') 
graph export "$output/TIFF/Lpoly_h_kw3.tif", as(tif) replace width(`res')
